###################################################################
## Replication code for "Gender in the Journals: A Dataset of Gender
## and Politics Research" by Barnett, FitzGerald, Krumbholz, and Lamba

# Code run using R version 4.0.3 (2020-10-10) and RStudio version 1.4.1103

library(Hmisc) # version 4.5-0
library(dplyr) # version 1.0.3
library(tidyr) # version 1.1.2
library(stringr) # version 1.4.0
library(ggplot2) # version 3.3.3

# Set working directory to local site of replication folder before proceeding
# setwd()

# Load summary data and restrict to 1980-2019 time frame
d_all <- readRDS("gender_research_summary.rds")
d <- filter(d_all, Year >= 1980)

# Set seed
set.seed(1234)



#### AGGREGATION/PROPORTION CALCULATION ####
## Journal Mean Proportions -----------------------------------------

## Sum journal-year counts, get mean non-duplicate keyword proportion,
# and use to estimate total unique articles by journal
agg <- d %>% 
    group_by(Journal, Abbreviation, Association, `Year Founded`,
             gender_journ, `Scopus H-Index`, `WOS H-Index`) %>% 
    summarise(across(count_sco:n_key_uniq, sum, na.rm = T),
              n_years_returned = sum(!is.na(prop_uniq_key)),
              n_years_total = n(),
              mean_prop_uniq_key = mean(prop_uniq_key, na.rm = T),
              LL95_prop_uniq_key = smean.cl.boot(prop_uniq_key)[[2]],
              UL95_prop_uniq_key = smean.cl.boot(prop_uniq_key)[[3]],
              est_total_journal = mean_prop_uniq_key * n_all,
              LL95_total_journal = LL95_prop_uniq_key * n_all,
              UL95_total_journal = UL95_prop_uniq_key * n_all)

## Apply to yearly data
d <- d %>% 
    left_join(agg[,c("Journal", "mean_prop_uniq_key",
                     "LL95_prop_uniq_key", "UL95_prop_uniq_key")]) %>% 
    mutate(est_total_journal = mean_prop_uniq_key * n_all,
           LL95_total_journal = LL95_prop_uniq_key * n_all,
           UL95_total_journal = UL95_prop_uniq_key * n_all)


## Pooled Mean Proportion -------------------------------------------

## Exclude journals with perfect duplication or non-duplication
outliers <- agg$Journal[agg$mean_prop_uniq_key %in% c(.5,1)]

## Get point estimate and 95% CI from all journal-years
pooled <- d %>% 
    filter(!(Journal %in% outliers)) %>% 
    group_by() %>% 
    summarise(mean = mean(prop_uniq_key, na.rm = T),
              LL95 = smean.cl.boot(prop_uniq_key)[[2]],
              UL95 = smean.cl.boot(prop_uniq_key)[[3]],
              n = sum(!is.na(prop_uniq_key)))
pooled # show pooled mean proportion estimate in console

## Apply to journal data aggregates
agg <- agg %>% 
    mutate(est_total_pooled = ifelse(!(Journal %in% outliers),
                                     pooled$mean * n_all,
                                     mean_prop_uniq_key * n_all),
           LL95_total_pooled = ifelse(!(Journal %in% outliers),
                                      pooled$LL95 * n_all,
                                      mean_prop_uniq_key * n_all),
           UL95_total_pooled = ifelse(!(Journal %in% outliers),
                                      pooled$UL95 * n_all,
                                      mean_prop_uniq_key * n_all))

## Apply to yearly data
d <- d %>% 
    mutate(est_total_pooled = ifelse(!(Journal %in% outliers),
                                     pooled$mean * n_all,
                                     mean_prop_uniq_key * n_all),
           LL95_total_pooled = ifelse(!(Journal %in% outliers),
                                      pooled$LL95 * n_all,
                                      mean_prop_uniq_key * n_all),
           UL95_total_pooled = ifelse(!(Journal %in% outliers),
                                      pooled$UL95 * n_all,
                                      mean_prop_uniq_key * n_all))



#### CODE CHOICE PROPORTIONS ####

## Using totals from pooled mean estimate ---------------------------

# Note that when dividing, the denominator for 95% CI is flipped
# since the upper/lower estimate of total articles will produce a
# smaller/larger estimate of the gender research proportion with a
# larger/smaller denominator.

agg <- agg %>% 
    mutate(`Unambig. Prop.` = Unambiguous / est_total_pooled,
           `Unambig. Prop. LL95` = Unambiguous / UL95_total_pooled,
           `Unambig. Prop. UL95` = Unambiguous / LL95_total_pooled,
           `Ambig. Prop.` = Ambiguous / est_total_pooled,
           `Ambig. Prop. LL95` = Ambiguous / UL95_total_pooled,
           `Ambig. Prop. UL95` = Ambiguous / LL95_total_pooled,
           `Combined Prop.` = `Unambig. Prop.` + `Ambig. Prop.`,
           `Combined Prop. LL95` = `Unambig. Prop. LL95` + `Ambig. Prop. LL95`,
           `Combined Prop. UL95` = `Unambig. Prop. UL95` + `Ambig. Prop. UL95`,
           `Ambig. Contribution` = Ambiguous / (Unambiguous + Ambiguous)) %>% 
    ungroup()


## Journal-Year Mean Proportions ------------------------------------

## Bayes estimator function
bayes <- function(n, y, mean, sd, max, c_int = .95) {
    
    ## Set hyperparameters and lambda grid
    alpha <- 1 + mean
    beta <- 1 + (1 / sd)
    lambda <- seq(0, max, .01)
    
    ## Get prior and posterior distribution
    prior <- dgamma(lambda, alpha, beta)
    posterior <- dgamma(lambda, alpha + y, beta + n)
    
    ## Numerical results
    bayes <- data.frame(lambda, prior, posterior)
    p_m <- (alpha + y) / (beta + n) # Posterior mean
    CI95 <- qgamma(c((1 - c_int) / 2, (1 + c_int) / 2),
                    alpha + y, beta + n) # 95% credible interval
    p_l <- bayes$lambda[bayes$posterior == 
                            max(bayes$posterior)] # Lambda at max prob
    
    return(list(p_m = p_m, LL95 = CI95[1], UL95 = CI95[2], p_l = p_l))
}

## Get prior info from summary data, pass to bayes estimator, and merge
agg <- d %>% 
    filter(gender_journ == 0) %>% 
    group_by(Abbreviation) %>% 
    summarise(n = n(),
              sum = sum(Combined),
              mean = mean(Combined),
              max = max(Combined)) %>% 
    rowwise() %>% 
    mutate(b_out = list(bayes(n, sum, mean(d$Combined[d$gender_journ == 0]),
                              sd(d$Combined[d$gender_journ == 0]), max))) %>% 
    unnest_wider(b_out) %>% 
    select(Abbreviation, p_m:p_l) %>% 
    left_join(agg, .) %>% 
    rename(`Combined Mean` = p_m,
           `Combined Mean LL95` = LL95,
           `Combined Mean UL95` = UL95,
           `Combined Prob. Mean` = p_l) %>% 
    mutate(`Comb. Mean Prop.` = `Combined Mean` / 
               (est_total_pooled / n_years_total),
           `Comb. Mean Prop LL95` = `Combined Mean LL95` / 
               (UL95_total_pooled / n_years_total),
           `Comb. Mean Prop UL95` = `Combined Mean UL95` / 
               (LL95_total_pooled / n_years_total))

## Average Yearly Counts --------------------------------------------
agg %>% # including PRQ outlier
    filter(gender_journ == 0) %>% 
    group_by() %>% 
    summarise(mean = round(mean(`Combined Mean`), 3),
              sd = round(sd(`Combined Mean`), 3)) # show in console
agg %>% # excluding PRQ outlier
    filter(gender_journ == 0, Abbreviation != "PRQ") %>% 
    group_by() %>% 
    summarise(mean = round(mean(`Combined Mean`), 3),
              sd = round(sd(`Combined Mean`), 3)) # show in console
agg %>% 
    filter(gender_journ == 0) %>% 
    select(Abbreviation, `Combined Mean`) %>% 
    arrange(-`Combined Mean`) %>% 
    as.data.frame() # show in console

## Combined prop. correlation with H-Index --------------------------
h_corr_table <- agg %>% 
    filter(gender_journ == 0) %>% 
    group_by() %>% 
    summarise(r_sco = cor.test(`Scopus H-Index`, `Combined Prop.`)[[4]][[1]],
              p_sco = cor.test(`Scopus H-Index`, `Combined Prop.`)[[3]],
              n_sco = sum(!is.na(`Scopus H-Index`)),
              r_wos = cor.test(`WOS H-Index`, `Combined Prop.`)[[4]][[1]],
              p_wos = cor.test(`WOS H-Index`, `Combined Prop.`)[[3]],
              n_wos = sum(!is.na(`WOS H-Index`)))
h_corr_table # show correlations in console



#### MAIN TABLES & FIGURES ####

## Main Table 1: Journals and Proportions in the Dataset ------------
main_tbl_1 <- agg %>% 
    filter(gender_journ == 0) %>% 
    select(Journal, Abbreviation, Combined, `Combined Prop.`) %>%
    mutate(`Combined Prop.` = round(`Combined Prop.`, 3)) %>% 
    arrange(Journal)
write.csv(main_tbl_1, "main_table_1.csv", row.names = F)


## Main Figure 1A: Overall Time Trend Graph --------------------------
main_fig_1a <- d %>% 
    group_by(gender_journ, Year) %>% 
    summarise(Combined = sum(Combined)) %>% 
    pivot_wider(names_from = gender_journ,
                names_prefix = "gj_",
                values_from = Combined) %>% 
    group_by(Year) %>% 
    summarise(`Non-gender-dedicated Only (n = 33)` = sum(gj_0, na.rm = T),
              `Includes Gender-dedicated (n = 37)` = sum(gj_1, na.rm = T)) %>% 
    mutate(`Includes Gender-dedicated (n = 37)` = 
               `Includes Gender-dedicated (n = 37)` +
               `Non-gender-dedicated Only (n = 33)`) %>% 
    pivot_longer(cols = c(`Non-gender-dedicated Only (n = 33)`,
                          `Includes Gender-dedicated (n = 37)`),
                 names_to = "Journal",
                 values_to = "Count") %>% 
    ggplot(aes(Year, Count, linetype = Journal)) +
    geom_line() +
    theme_bw() +
    theme(legend.position = c(.3, .8),
          legend.background = element_rect(color = "black"))

## Save plot
ggsave(file.path(getwd(), "main_fig_1a.tiff"),
       main_fig_1a,
       width = 5.5,
       height = 3.5,
       units = "in",
       dpi = 1200,
       compression = "lzw")

## Main Figure 1B: Time Trend Graph with Proportions --------------------------

# All journals
all_props <- d %>% 
    group_by(Year) %>% 
    summarise(Combined = sum(Combined),
              Est_pooled = sum(est_total_pooled)) %>% 
    mutate(Proportion = Combined / Est_pooled,
           Journal = "Includes Gender-dedicated (n = 37)")

# Non-gender only
ng_props <- d %>% 
    filter(gender_journ == 0) %>%
    group_by(Year) %>% 
    summarise(Combined = sum(Combined),
              Est_pooled = sum(est_total_pooled)) %>% 
    mutate(Proportion = Combined / Est_pooled,
           Journal = "Non-gender-dedicated Only (n = 33)") 

# Combine data and plot
main_fig_1b <- bind_rows(all_props, ng_props) %>%
    ggplot(aes(Year, Proportion, linetype = Journal)) +
        geom_line() +
        scale_y_continuous(limits = c(0, 0.25)) +
        theme_bw() +
        theme(legend.position = c(.3, .8),
              legend.background = element_rect(color = "black")) 

## Save plot
ggsave(file.path(getwd(), "main_fig_1b.tiff"),
       main_fig_1b,
       width = 5.5,
       height = 3.5,
       units = "in",
       dpi = 1200,
       compression = "lzw")



## Main Figure 2A: Time trends for "recent uptick" journals -----------
main_fig_2a <- d %>% 
    filter(Abbreviation %in% c("AJPS","BJPS","CPS","JOP",
                               "NPS","POP","PB","PRQ")) %>% 
    ggplot(aes(Year, Combined)) +
    geom_line() +
    scale_y_continuous(limits = c(0, 20)) +
    labs(y = "Count") +
    facet_wrap(~ Journal, ncol = 2) +
    theme_bw()

## Save plot
ggsave(file.path(getwd(), "main_fig_2a.tiff"),
       main_fig_2a,
       width = 6.5,
       height = 4.5,
       units = "in",
       dpi = 1200,
       compression = "lzw")

## Main Figure 2B: Time trends for "recent uptick" journals using props -----------
main_fig_2b <- d %>% 
    filter(Abbreviation %in% c("AJPS","BJPS","CPS","JOP",
                               "NPS","POP","PB","PRQ")) %>% 
    group_by(Journal, Year) %>%
    summarise(Combined = sum(Combined),
              Est_pooled = sum(est_total_pooled)) %>% 
    mutate(Proportion = Combined / Est_pooled) %>%
    ggplot(aes(Year, Proportion)) +
    geom_line() +
    scale_y_continuous(limits = c(0, 0.4)) +
    labs(y = "Proportion") +
    facet_wrap(~ Journal, ncol = 2) +
    theme_bw()

## Save plot
ggsave(file.path(getwd(), "main_fig_2b.tiff"),
       main_fig_2b,
       width = 6.5,
       height = 4.5,
       units = "in",
       dpi = 1200,
       compression = "lzw")


#### APPENDIX TABLES & FIGURES ####

## Appendix Table A1: Journals in the Dataset -----------------------
append_tbl_A1 <- agg %>% 
    select(Journal:`Year Founded`, `Scopus H-Index`, `WOS H-Index`) %>%
    rename(Founded = `Year Founded`,
           `SCOPUS H-I` = `Scopus H-Index`,
           `WOS H-I` = `WOS H-Index`) %>% 
    arrange(Journal)
write.csv(append_tbl_A1, "append_table_A1.csv", row.names = F)


## Appendix Table A2: Coding Decisions by Journal -------------------
append_tbl_A2 <- d_all %>% 
    group_by(Abbreviation) %>% 
    summarise(across(c(Unambiguous, Ambiguous,
                       `Exclude (Content)`, `Exclude (Type)`),
                     sum, na.rm = T)) %>% 
    mutate(Total = Unambiguous + Ambiguous +
               `Exclude (Content)` + `Exclude (Type)`) %>% 
    rename(Journal = Abbreviation)
write.csv(append_tbl_A2, "append_table_A2.csv", row.names = F)


## Appendix Table A3: Search String Validation ----------------------
append_tbl_A3 <- tibble(Journal = c("IFJP", "JWPP", "P&G", "SP", "W&P"),
                        sco_keywords = c(301, 229, 295, 327, 440),
                        sco_non_key = c(333, 230, 301, 348, 481),
                        sco_pct = round(sco_keywords / sco_non_key, 3),
                        wos_keywords = c(350, 263, 314, 444, 260),
                        wos_non_key = c(393, 264, 320, 475, 267),
                        wos_pct = round(wos_keywords / wos_non_key, 3))
write.csv(append_tbl_A3, "append_table_A3.csv", row.names = F)


## Appendix Figures E1A-3A: Time Trends for Additional Journals -------

## Include only non-gender journals with >= 10 years of data and that are
# not already included in the main manuscript
g_journals <- d %>% 
    filter(gender_journ == 0,
           !(Abbreviation %in% c("AJPS","BJPS","CPS","JOP",
                                 "NPS","POP","PB","PRQ")),
           Journal %in% agg$Journal[agg$n_years_returned >= 10]) %>% 
    arrange(Journal) %>% 
    mutate(Journal = ifelse(Journal == "JOURNAL OF HEALTH POLITICS, POLICY AND LAW",
                            "JHPPL",
                            ifelse(Journal == "JOURNAL OF INFORMATION TECHNOLOGY AND POLITICS",
                                   "JITP",
                                   Journal))
           )

## Plot in 3 sets of 6
append_fig_E1a <-g_journals %>% 
    filter(Journal %in% unique(Journal)[1:6]) %>% 
    ggplot(aes(Year, Combined)) +
    geom_line() +
    scale_y_continuous(limits = c(0, 20)) +
    labs(y = "Count") +
    facet_wrap(~ Journal, ncol = 2) +
    theme_bw()

append_fig_E2a <-g_journals %>% 
    filter(Journal %in% unique(Journal)[7:12]) %>% 
    ggplot(aes(Year, Combined)) +
    geom_line() +
    scale_y_continuous(limits = c(0, 20)) +
    labs(y = "Count") +
    facet_wrap(~ Journal, ncol = 2) +
    theme_bw()

append_fig_E3a <-g_journals %>% 
    filter(Journal %in% unique(Journal)[13:18]) %>% 
    ggplot(aes(Year, Combined)) +
    geom_line() +
    scale_y_continuous(limits = c(0, 20)) +
    labs(y = "Count") +
    facet_wrap(~ Journal, ncol = 2) +
    theme_bw()

## Save plots
append_fig_Ea <- list(append_fig_E1a = append_fig_E1a,
                     append_fig_E2a = append_fig_E2a,
                     append_fig_E3a = append_fig_E3a)
lapply(seq_along(append_fig_Ea),
       function(i) {
           ggsave(file.path(getwd(), paste0(names(append_fig_Ea)[[i]], ".tiff")),
                  append_fig_Ea[[i]],
                  width = 6.5,
                  height = 4.5,
                  units = "in",
                  dpi = 1200,
                  compression = "lzw")
       })


## Appendix Figures E1B-3B: Time Trends for Additional Journals (Proportions) -------

## Include only non-gender journals with >= 10 years of data and that are
# not already included in the main manuscript
g_journals <- d %>% 
    filter(gender_journ == 0,
           !(Abbreviation %in% c("AJPS","BJPS","CPS","JOP",
                                 "NPS","POP","PB","PRQ")),
           Journal %in% agg$Journal[agg$n_years_returned >= 10]) %>% 
    arrange(Journal) %>% 
    mutate(Journal = ifelse(Journal == "JOURNAL OF HEALTH POLITICS, POLICY AND LAW",
                            "JHPPL",
                            ifelse(Journal == "JOURNAL OF INFORMATION TECHNOLOGY AND POLITICS",
                                   "JITP",
                                   Journal)),
           Proportion = Combined / est_total_pooled
    )

## Plot in 3 sets of 6
append_fig_E1b <-g_journals %>% 
    filter(Journal %in% unique(Journal)[1:6]) %>% 
    ggplot(aes(Year, Proportion)) +
    geom_line() +
    scale_y_continuous(limits = c(0, 0.4)) +
    labs(y = "Proportion") +
    facet_wrap(~ Journal, ncol = 2) +
    theme_bw()

append_fig_E2b <-g_journals %>% 
    filter(Journal %in% unique(Journal)[7:12]) %>% 
    ggplot(aes(Year, Proportion)) +
    geom_line() +
    scale_y_continuous(limits = c(0, 0.4)) +
    labs(y = "Proportion") +
    facet_wrap(~ Journal, ncol = 2) +
    theme_bw()

append_fig_E3b <-g_journals %>% 
    filter(Journal %in% unique(Journal)[13:18]) %>% 
    ggplot(aes(Year, Proportion)) +
    geom_line() +
    scale_y_continuous(limits = c(0, 0.7)) + # NOTE different y-axis
    labs(y = "Proportion") +
    facet_wrap(~ Journal, ncol = 2) +
    theme_bw()

## Save plots
append_fig_Eb <- list(append_fig_E1b = append_fig_E1b,
                     append_fig_E2b = append_fig_E2b,
                     append_fig_E3b = append_fig_E3b)
lapply(seq_along(append_fig_Eb),
       function(i) {
           ggsave(file.path(getwd(), paste0(names(append_fig_Eb)[[i]], ".tiff")),
                  append_fig_Eb[[i]],
                  width = 6.5,
                  height = 4.5,
                  units = "in",
                  dpi = 1200,
                  compression = "lzw")
       })



## Appendix Figure E4A: Gender-Dedicated Journals Time Trends --------
append_fig_E4a <- d %>% 
    filter(gender_journ == 1) %>% 
    mutate(Journal = ifelse(Journal %in% c("WOMEN AND POLITICS",
                                           "JOURNAL OF WOMEN, POLITICS AND POLICY"),
                            "WOMEN AND POLITICS / JWPP",
                            ifelse(Journal == "INTERNATIONAL FEMINIST JOURNAL OF POLITICS",
                                   "IFJP", Journal))) %>% 
    ggplot(aes(Year, Combined)) +
    stat_summary(fun = sum, geom = "line") +
    scale_y_continuous(limits = c(0, 60)) +
    labs(y = "Count") +
    facet_wrap(~ Journal, ncol = 2) +
    theme_bw()

## Save plot
ggsave(file.path(getwd(), "append_fig_E4a.tiff"),
       append_fig_E4a,
       width = 6.5,
       height = 4.5,
       units = "in",
       dpi = 1200,
       compression = "lzw")

## Appendix Figure E4B: Gender-Dedicated Journals Time Trends (Proportions) --------
append_fig_E4b <- d %>% 
    filter(gender_journ == 1) %>% 
    mutate(Journal = ifelse(Journal %in% c("WOMEN AND POLITICS",
                                           "JOURNAL OF WOMEN, POLITICS AND POLICY"),
                            "WOMEN AND POLITICS / JWPP",
                            ifelse(Journal == "INTERNATIONAL FEMINIST JOURNAL OF POLITICS",
                                   "IFJP", Journal)),
           Proportion = Combined / est_total_pooled) %>% 
    ggplot(aes(Year, Proportion)) +
    stat_summary(fun = sum, geom = "line") +
    scale_y_continuous(limits = c(0, 1)) +
    labs(y = "Proportion") +
    facet_wrap(~ Journal, ncol = 2) +
    theme_bw()

## Save plot
ggsave(file.path(getwd(), "append_fig_E4b.tiff"),
       append_fig_E4b,
       width = 6.5,
       height = 4.5,
       units = "in",
       dpi = 1200,
       compression = "lzw")



## Appendix Figure E5: Average Yearly Proportions -------------------

## Get pooled mean and 95% CI
pooled <- agg %>% 
    filter(gender_journ == 0) %>% 
    group_by() %>% 
    summarise(n = n(),
              mean = round(mean(`Comb. Mean Prop.`), 3),
              LL95 = round(mean(t.test(`Comb. Mean Prop.`)[[4]][[1]]), 3),
              UL95 = round(mean(t.test(`Comb. Mean Prop.`)[[4]][[2]]), 3))

## Order journals by mean, descending
f_order <- agg %>% 
    arrange(`Comb. Mean Prop.`) %>% 
    mutate(lbl = paste0(Abbreviation, " (", n_years_total, ")")) %>% 
    select(Abbreviation, lbl)

## Plot
append_fig_E5 <- agg %>% 
    filter(gender_journ == 0) %>% 
    mutate(Abbreviation = factor(Abbreviation,
                                 levels = f_order$Abbreviation,
                                 labels = f_order$lbl)) %>% 
    ggplot(aes(Abbreviation, `Comb. Mean Prop.`)) +
    geom_point(stat = "identity", size = 2) +
    geom_errorbar(aes(ymin = `Comb. Mean Prop LL95`,
                      ymax = `Comb. Mean Prop UL95`)) +
    geom_hline(yintercept = pooled$mean, linetype = "solid") +
    geom_hline(yintercept = c(pooled$LL95, pooled$UL95), linetype = "dashed") +
    labs(x = "Journal (n Years)", y = "Mean Proportion (95% CI)") +
    coord_flip() +
    theme_bw()

## Save plot
ggsave(file.path(getwd(), "append_fig_E5.tiff"),
       append_fig_E5,
       width = 5.5,
       height = 4.5,
       units = "in",
       dpi = 1200,
       compression = "lzw")


## Appendix Figure E6: Average Yearly Counts ------------------------

## Get pooled mean and 95% CI
pooled <- agg %>% 
    filter(gender_journ == 0) %>% 
    group_by() %>% 
    summarise(n = n(),
              mean = round(mean(`Combined Mean`), 3),
              LL95 = round(t.test(`Combined Mean`)[[4]][[1]], 3),
              UL95 = round(t.test(`Combined Mean`)[[4]][[2]], 3))

## Order journals by mean, descending
f_order <- agg %>% 
    arrange(`Combined Mean`) %>% 
    mutate(lbl = paste0(Abbreviation, " (", n_years_total, ")")) %>% 
    select(Abbreviation, lbl)

## Plot
append_fig_E6 <- agg %>% 
    filter(gender_journ == 0) %>% 
    mutate(Abbreviation = factor(Abbreviation,
                                 levels = f_order$Abbreviation,
                                 labels = f_order$lbl)) %>% 
    ggplot(aes(Abbreviation, `Combined Mean`)) +
    geom_point(stat = "identity", size = 2) +
    geom_errorbar(aes(ymin = `Combined Mean LL95`,
                      ymax = `Combined Mean UL95`)) +
    geom_hline(yintercept = pooled$mean, linetype = "solid") +
    geom_hline(yintercept = c(pooled$LL95, pooled$UL95), linetype = "dashed") +
    labs(x = "Journal (n Years)", y = "Mean Count (95% CI)") +
    coord_flip() +
    theme_bw()

## Save plot
ggsave(file.path(getwd(), "append_fig_E6.tiff"),
       append_fig_E6,
       width = 5.5,
       height = 4.5,
       units = "in",
       dpi = 1200,
       compression = "lzw")



#end